#include "MModel.h"

#include <Eigen/Eigenvalues>

#include "Utility/Algebra/ODESystem.h"
#include "Utility/Algebra/sgArpack.h"
#include "Utility/Algebra/sgPardiso.h"
#include "Utility/Logging/SGLogger.h"
#ifdef SG_LICENSE_CHECK
#include <Utility/Security/LicenseCheck.h>
#endif

#include "DataStructure/Input/Constraint/SPCADDData.h"
#include "DataStructure/Input/Dynamic/InitConditionType.h"
#include "DataStructure/Input/Parameter/TSTEPData.h"
#include "DataStructure/Input/Solver/DOFDirections.h"

#include "DBManager/IService/Input/Constraint/IConstraintSummaryService.h"
#include "DBManager/IService/Input/Constraint/ISPCADDService.h"
#include "DBManager/IService/Input/Coordinate/ICORDSummaryService.h"
#include "DBManager/IService/Input/Coordinate/ITRANSFORMATIONService.h"
#include "DBManager/IService/Input/Dynamic/IFREQ1Service.h"
#include "DBManager/IService/Input/Dynamic/ITICService.h"
#include "DBManager/IService/Input/Element/IElementService.h"
#include "DBManager/IService/Input/Load/IDLOADService.h"
#include "DBManager/IService/Input/Load/ILoadSummaryService.h"
#include "DBManager/IService/Input/Load/ILoadCombinationService.h"

#include "DBManager/IService/Input/Load/IRLOAD1Service.h"
#include "DBManager/IService/Input/Load/IRLOAD2Service.h"
#include "DBManager/IService/Input/Load/ITLOAD1Service.h"
#include "DBManager/IService/Input/Load/ITLOAD2Service.h"
#include "DBManager/IService/Input/Parameter/IParameterService.h"
#include "DBManager/IService/Input/Parameter/ITSTEPService.h"
#include "DBManager/IService/Input/Subcase/ISubcaseService.h"
#include "DBManager/IService/Input/Table/ITABDMP1Service.h"
#include "DBManager/IService/Result/Nodal/IACCELERATIONService.h"
#include "DBManager/IService/Result/Nodal/IACCELERATION_CPLXService.h"
#include "DBManager/IService/Result/Nodal/IDISPLACEMENTService.h"
#include "DBManager/IService/Result/Nodal/IDISPLACEMENT_CPLXService.h"
#include "DBManager/IService/Result/Nodal/IEIGENVECTORService.h"
#include "DBManager/IService/Result/Nodal/IEIGENVECTOR_CPLXService.h"
// #include "DBManager/IService/Result/IResultDomainsService.h"
#include "DBManager/IService/Result/Nodal/IVELOCITYService.h"
#include "DBManager/IService/Result/Nodal/IVELOCITY_CPLXService.h"
#include "DBManager/IService/Result/Summary/IEIGENVALUEService.h"
#include "DBManager/IService/Result/Summary/IEIGENVALUE_CPLXService.h"
#include "DBManager/IService/Result/Random/ICROSS_PSDService.h"

#include "SGFem/Pre/Assembly/CoordinateProcessor.h"
#include "SGFem/Pre/Assembly/DataAdapter.h"

#include "SGFem/ElementCalculator/Factory/ElementCalculatorFactory.h"
#include "SGFem/Task/Common/LoadCalculator.h"

#include "SGFem/Constraint/RElement/BoundaryCondition.h"
// #include "DBManager/IService/Input/Element/ICHBDYPService.h"
// #include "DBManager/IService/Input/Element/ICHBDYGService.h"
#include "SGFem/ElementCalculator/ThermalElement/CHBDYGAREACalculator.h"
#include "SGFem/ElementCalculator/ThermalElement/CHBDYPLINECalculator.h"

using SG::Algebra::Complex_t;
using namespace SG::DataStructure::Common;
using namespace SG::DataStructure::FEM;
using namespace SG::DataStructure::Post;
using namespace SG::DBManager;

using SG::FEMSolver::EigenVector;
using SG::FEMSolver::AssemblyEigen;
using SG::FEMSolver::Constraint::BoundaryCondition;
using SG::FEMSolver::Element::ElementCalculatorFactory;

using namespace Kernel;

MModel::MModel (SG::DataStructure::FEM::SOLType analysisType, SG::DBManager::DBServiceFactorySP dbServiceFactory)
    : m_dbServiceFactory (dbServiceFactory), m_analysisType (analysisType)
{}

MModel::~MModel () = default;

void MModel::TimeDomainLoadInit (const SG::DataStructure::FEM::DLOADData& dloadData, SG::DataStructure::Common::Id loadSetId)
{
    m_timeDynaLoadSet.clear ();
    for (const auto& sl : dloadData.m_SLList)
    {
        // Query TLOAD1
        std::vector<SG::DataStructure::FEM::TLOAD1Data> tload1List;
        if (Status::STATUS_Success == m_dbServiceFactory->get<ITLOAD1Service> ()->FindBySId (sl.m_LI, tload1List))
        {
            for (const auto& tload1 : tload1List)
            {
                // 查找单元静力载荷集编号
                auto loadIds = getStaticLoadId (tload1, m_dbServiceFactory, loadSetId);
                for (auto sid : loadIds)
                {
                    m_loadCases.insert (sid);
                }

                auto calc = std::make_shared<TLoad1Calc> (tload1, loadSetId, sl.m_SI * dloadData.m_S);
                m_timeDynaLoadSet.push_back (calc);
            }
        }

        // Query TLOAD2
        std::vector<SG::DataStructure::FEM::TLOAD2Data> tload2List;
        if (Status::STATUS_Success == m_dbServiceFactory->get<ITLOAD2Service> ()->FindBySId (sl.m_LI, tload2List))
        {
            for (const auto& tload2 : tload2List)
            {
                auto loadIds = getStaticLoadId (tload2, m_dbServiceFactory, loadSetId);
                for (auto sid : loadIds)
                {
                    m_loadCases.insert (sid);
                }
                auto calc = std::make_shared<TLoad2Calc> (tload2, loadSetId, sl.m_SI * dloadData.m_S);
                m_timeDynaLoadSet.push_back (calc);
            }
        }
    }
}

void MModel::FreqDomainLoadInit (const SG::DataStructure::FEM::DLOADData& dloadData, SG::DataStructure::Common::Id loadSetId)
{
    m_freqDynaLoadSet.clear ();
    for (const auto& sl : dloadData.m_SLList)
    {
        // Query RLOAD1
        std::vector<SG::DataStructure::FEM::RLOAD1Data> rload1List;
        if (Status::STATUS_Success == m_dbServiceFactory->get<IRLOAD1Service> ()->FindBySId (sl.m_LI, rload1List))
        {
            for (const auto& rload1 : rload1List)
            {
                auto loadIds = getStaticLoadId (rload1, m_dbServiceFactory, loadSetId);
                for (auto sid : loadIds)
                {
                    m_loadCases.insert (sid);
                }
                auto calc = std::make_shared<RLoad1Calc> (rload1, loadSetId, sl.m_SI * dloadData.m_S);
                m_freqDynaLoadSet.push_back (calc);
            }
        }

        // Query RLOAD2
        std::vector<SG::DataStructure::FEM::RLOAD2Data> rload2List;
        if (Status::STATUS_Success == m_dbServiceFactory->get<IRLOAD2Service> ()->FindBySId (sl.m_LI, rload2List))
        {
            for (const auto& rload2 : rload2List)
            {
                auto loadIds = getStaticLoadId (rload2, m_dbServiceFactory, loadSetId);
                for (auto sid : loadIds)
                {
                    m_loadCases.insert (sid);
                }
                auto calc = std::make_shared<RLoad2Calc> (rload2, loadSetId, sl.m_SI * dloadData.m_S);
                m_freqDynaLoadSet.push_back (calc);
            }
        }

        // TODO(29794): 频响可能使用 TLOAD1， TLOAD2载荷。 需要进行扩展
    }
}

void MModel::Initialize (const SG::DataStructure::Common::IdArray& subcaseIds)
{
// 检查授权
#ifdef SG_LICENSE_CHECK
    SG::Utility::SGCheckLicense ();
#endif

    // subcase 为空
    if (0 == subcaseIds.size ())
    {
        return;
    }
    auto& assembly = SG::FEMSolver::AssemblyEigen::GetInstance ();
    assembly.Clear ();

    // 获取Subcase
    // 所有工况对应的边界条件都相同， 因此只选择第一个工况
    auto                            para = m_dbServiceFactory->get<ISubcaseService> ()->FindSubCase (subcaseIds[0]);
    SG::DataStructure::FEM::Subcase subcase (para);

    // 获取静载荷信息
    // std::set<SG::DataStructure::Common::Id> loadCases;
    if (SOLType::SOL_SESTATIC == m_analysisType || SOLType::SOL_NLSCSH == m_analysisType)
    {
        m_loadCases = m_dbServiceFactory->get<ISubcaseService> ()->GetLoadsetIds (subcaseIds);
        // 热分析把边界条件视为载荷热当成载荷[同一组工况，边界条件相同]，因此只需要处理一次
        if (SOLType::SOL_NLSCSH == m_analysisType)
        {
            auto spcId = subcase.GetSPCId ();
            if (0 != spcId)
            {
                m_loadCases.emplace (spcId);
            }
        }
    }
    else
    {
        // 获取 DLOAD 数据
        SG::DataStructure::FEM::DLOADData dloadData;
        auto status = m_dbServiceFactory->get<IDLOADService> ()->FindBySId (subcase.GetDLoadId (), _OUT dloadData);
        // 如果查不到，就构造一个默认 DLOAD
        if (status == SG::DataStructure::Common::Status::STATUS_Failure)
        {
            dloadData.m_S = 1.0;
            dloadData.m_SLList.push_back ({ 1.0, subcase.GetDLoadId () });
        }

        auto loadSetId = subcase.GetLoadsetId ();

        // 遍历 DLOAD 中的 SLList，查找对应的 动力载荷
        // 并且查询用到的静力载荷id
        if (SOLType::SOL_SEDTRAN == m_analysisType || SOLType::SOL_SEMTRAN == m_analysisType)
        {
            TimeDomainLoadInit (dloadData, loadSetId);
        }

        if (SOLType::SOL_SEDFREQ == m_analysisType || SOLType::SOL_SEMFREQ == m_analysisType)
        {
            FreqDomainLoadInit (dloadData, loadSetId);
        }
    }

    // 完成节点和自由度排序。开辟整体矩阵及向量的存储空间
    assembly.Initialize (m_dbServiceFactory, subcase, m_loadCases);

    // 激活质量矩阵计算 (非静力，非热力学工况 或者具有重力载荷)
    const auto GRAVDatas = m_dbServiceFactory->get<ILoadSummaryService> ()->FindGRAV (m_loadCases);
    if ((SOLType::SOL_SESTATIC != m_analysisType) || GRAVDatas.size () > 0)
    {
        assembly.ActiveMassMatrix ();
    }

    // 初始化，计算所有局部坐标系的基矢量和全局坐标系
    auto& coordHandler = SG::FEMSolver::CoordinateProcessor::GetInstance ();
    coordHandler.Initialize (m_dbServiceFactory);

    // 存储TRANSFORMATION
    std::vector<SG::DataStructure::Common::Id> allCoordId   = m_dbServiceFactory->get<ICORDSummaryService> ()->GetAllId ();
    auto                                       transService = m_dbServiceFactory->get<ITRANSFORMATIONService> ();

    std::vector<SG::DataStructure::FEM::TRANSFORMATIONData> transDataList;
    SG::DataStructure::FEM::TRANSFORMATIONData              transData;
    for (auto& cid : allCoordId)
    {
        transData.m_RDATAList.clear ();
        transData.m_RDATAList.resize (12);  // 12个实数数据

        const SG::DataStructure::FEM::ICoordinateP coord = coordHandler.GetCoord (cid);

        transData.m_id = coord->m_id;

        transData.setType (coord->m_type);

        transData.setOrigin (coord->GetOrigin ());
        transData.setBaseX (coord->GetBaseX ());
        transData.setBaseY (coord->GetBaseY ());
        transData.setBaseZ (coord->GetBaseZ ());

        transDataList.emplace_back (transData);
    }
    if (transDataList.size () > 0)
    {
        // 将坐标系数据存储到数据库中
        if (transService->MAdd (transDataList) != Status::STATUS_Success)
        {
            std::cerr << "Failed to add transformation data to the database.\n";
            return;
        }
    }
}

// 静力分析，计算每个单元的刚度阵和质量阵并组装
Status MModel::CalcMatrixEquation (int subcaseId)
{
    std::cout << "===== compute stiffness \n";
    const auto elementTypes = m_dbServiceFactory->get<IElementService> ()->GetAllElementTypes ();
    for (const auto& type : elementTypes)
    {
        if ((type == SG::DataStructure::FEM::ELETYPE_CHBDYP) || type == SG::DataStructure::FEM::ELETYPE_CHBDYG)
        {
            continue;
        }

        auto eleCalculator = ElementCalculatorFactory::Create (type);
        if (NULL == eleCalculator)
        {
            std::cout << "can 't find element calc " << std::endl;
        }

        eleCalculator->Initialize (type);
        eleCalculator->ComputeElementMatries (m_dbServiceFactory);
    }
    std::cout << "===== compute stiffness finish\n";
    auto&   assembly = AssemblyEigen::GetInstance ();
    auto    para     = m_dbServiceFactory->get<ISubcaseService> ()->FindSubCase (subcaseId);
    Subcase subcase (para);

    if (SOLType::SOL_NLSCSH == m_analysisType)
    {
        // 计算热边界，对流
        auto         spcId = subcase.GetSPCId ();
        std::set<Id> spcIdSet;
        const auto   spcTypes = m_dbServiceFactory->get<IConstraintSummaryService> ()->GetType (spcId);

        if (0 != std::count (spcTypes.begin (), spcTypes.end (), ConstraintType::CONTYPE_SPCADD))
        {
            std::vector<SPCADDData> spcAdddatas;
            m_dbServiceFactory->get<ISPCADDService> ()->FindBySId (spcId, spcAdddatas);
            for (const auto& spcAdd : spcAdddatas)
            {
                for (const auto& spdIdJ : spcAdd.m_subSets)
                {
                    spcIdSet.insert (spdIdJ);
                }
            }
        }
        else
        {
            spcIdSet.insert (spcId);
        }

        for (const auto& id : spcIdSet)
        {
            SG::FEMSolver::Element::CHBDYPLINECalculator bdyeleCalculator;
            bdyeleCalculator.Initialize (SG::DataStructure::FEM::ELETYPE_CHBDYP);
            bdyeleCalculator.ComputeConvMatries (m_dbServiceFactory, id);

            SG::FEMSolver::Element::CHBDYGAREACalculator bdyeleCalculator2;
            bdyeleCalculator2.Initialize (SG::DataStructure::FEM::ELETYPE_CHBDYG);
            bdyeleCalculator2.ComputeConvMatries (m_dbServiceFactory, id);
        }
    }

    assembly.Assembly (m_dbServiceFactory);

    // 计算整体载荷列向量
    LoadCalculator::ComputeStaticLoad (m_dbServiceFactory, m_loadCases);

    ///< 计算mpc方程对应的控制矩阵
    BoundaryCondition::ComputeControlMatrix (m_dbServiceFactory, subcase);

    assembly.Constraint (m_dbServiceFactory, subcaseId);
    std::cout << "== assembly finish" << std::endl;

    return Status::STATUS_Success;
}

void MModel::ComputeLoadVector(const Id loadId, EigenVector& f) const
{
    auto& assembly = AssemblyEigen::GetInstance();
    LoadCombination load;
    const auto isLoadCase = m_dbServiceFactory -> get<ILoadCombinationService>() -> FindById (loadId, load);
    if (Status::STATUS_Success == isLoadCase)
    {
        //SG::DataStructure::FEM::print (load);
        assembly.ComputeBVector (m_dbServiceFactory, load, _OUT f);
    }
    else
    {
        assembly.ComputeBVector (loadId,  _OUT f);
    }
}

void MModel::ComputeStiffnessMatrix (SG::Algebra::SpMat_t& K)
{
    K = SG::FEMSolver::AssemblyEigen::GetInstance ().GetStiffnessMatrix ();
}

void MModel::ComputeMassMatrix (SG::Algebra::SpMat_t& M)
{
    M                              = SG::FEMSolver::AssemblyEigen::GetInstance ().GetCoupledMassMatrix ();
    SG::Utility::SGParameter param = m_dbServiceFactory->get<IParameterService> ()->GetSolverParameter ();
    if (param.Has ("WTMASS"))
    {
        auto coeff = std::stod (param["WTMASS"].GetString ());
        M          = M * coeff;
    }
}

void MModel::ComputeLumpedMassMatrix (SG::Algebra::Vec_t& M)
{
    SGWarn << __FILE__ << ":" << __LINE__ << " NOT IMPLEMENTED";
}

void MModel::ComputeDampingMatrix (SG::Algebra::SpMat_t& C)
{
    SG::Utility::SGParameter para = m_dbServiceFactory->get<IParameterService> ()->GetSolverParameter ();

    auto&       assembly = SG::FEMSolver::AssemblyEigen::GetInstance ();
    std::size_t n        = assembly.GetFreeDOFSize ();
    C                    = SG::Algebra::SpMat_t (n, n);
    C.setZero ();

    // 瑞利阻尼
    if (para.Has ("ALPHA1"))
    {
        Real alpha{ 0.0 };
        alpha         = std::stod (para["ALPHA1"]["V1"].GetString ());
        const auto& M = SG::FEMSolver::AssemblyEigen::GetInstance ().GetCoupledMassMatrix ();

        C += alpha * M;
    }

    if (para.Has ("ALPHA2"))
    {
        Real beta{ 0.0 };
        beta          = std::stod (para["ALPHA2"]["V1"].GetString ());
        const auto& K = SG::FEMSolver::AssemblyEigen::GetInstance ().GetStiffnessMatrix ();
        C += beta * K;
    }

    // 添加结构阻尼
    // 如果同时存在G和W3参数
    if (para.Has ("G") && para.Has ("W3"))
    {
        Real DampG         = 0.0;
        Real DampW3        = 1E10;
        DampG              = std::stod (para["G"].GetString ());
        DampW3             = std::stod (para["W3"].GetString ());
        const auto& K      = SG::FEMSolver::AssemblyEigen::GetInstance ().GetStiffnessMatrix ();
        const Real  factor = DampG / DampW3;
        C += factor * K;
    }

    // ToDo添加其他的阻尼
    // C += 阻尼单元 + B2GG +B2PP;

    // ToDo单元结构阻尼->等效粘性阻尼
}

SG::Algebra::Real_t MModel::ComputeDampingRatio (SG::Algebra::Real_t omega)
{
    auto                            para = m_dbServiceFactory->get<ISubcaseService> ()->FindSubCase (1);
    SG::DataStructure::FEM::Subcase subcase (para);
    const auto                      sdampId = subcase.GetSdampingId ();

    // 如果没有设置结构阻尼，直接返回Xi = 0
    if (sdampId == 0)
    {
        return 0.0;
    }

    SG::DataStructure::FEM::TABDMP1Data TABDMP1Data;
    auto                                status = m_dbServiceFactory->get<ITABDMP1Service> ()->FindById (sdampId, TABDMP1Data);

    // 如果获取失败或数据为空，返回Xi = 0
    if (Status::STATUS_Failure == status || TABDMP1Data.m_DAMPList.empty ())
    {
        return 0.0;
    }

    // 如果成功获取到数据，进行插值计算，TABDMP1输出自动转换为阻尼系数G
    return TABDMP1Data.Interpolate (omega) * 0.5;
}

// void MModel::ComputeForce (SG::Algebra::Vec_t& F)
// {
//     SGWarn << __FILE__ << ":" << __LINE__ << " NOT IMPLEMENTED";
// }

Status MModel::QueryTimeStep (Id tStepId, std::vector<Real>& tCompute, std::vector<Real>& tOut) const
{
    // 计算时间步
    SG::DataStructure::FEM::TSTEPData TSTEPData;
    if (Status::STATUS_Success == m_dbServiceFactory->get<ITSTEPService> ()->FindBySId (tStepId, TSTEPData))
    {
        tCompute = TSTEPData.GetComputeTimeVector ();
    }
    else
    {
        SGError << "can't find TSTEPData setting";
        return Status::STATUS_Failure;
    }

    return Status::STATUS_Success;
}

void MModel::ComputeInitialCondition (Id                                       icId,
                                      SG::DataStructure::FEM::InitConditonType type,
                                      SG::Algebra::Vec_t&                      U,
                                      SG::Algebra::Vec_t&                      V) const
{
    // 初始条件默认为0
    auto&       assembly = SG::FEMSolver::AssemblyEigen::GetInstance ();
    std::size_t nDof     = assembly.GetFreeDOFSize ();
    U                    = SG::Algebra::Vec_t (nDof);
    U.setZero ();
    V = SG::Algebra::Vec_t (nDof);
    V.setZero ();

    std::vector<SG::DataStructure::FEM::TICData> TICData;
    if (Status::STATUS_Failure == m_dbServiceFactory->get<ITICService> ()->FindBySId (icId, TICData))
    {
        return;
    }

    if (InitConditonType::INITCONDTYPE_PHYSICAL == type)
    {
        std::size_t nodeIndex;
        int         dofGlobalIds;

        for (const auto& tic : TICData)
        {
            assembly.FindNodeIndex (tic.m_G, _OUT nodeIndex);
            assembly.GetGlobalDOFId (nodeIndex, tic.GetDirection (), _OUT dofGlobalIds);
            U[dofGlobalIds] = tic.m_U0;
            V[dofGlobalIds] = tic.m_V0;
        }
    }
    else
    {
        SGError << "can't support this init condition ";
    }
}

void MModel::ComputeDerivatives (const SG::Algebra::Real_t& t,
                                 const SG::Algebra::Vec_t&  u,
                                 const SG::Algebra::Vec_t&  v,
                                 SG::Algebra::Vec_t& _OUT   a)
{
    SG::Algebra::SpMat_t K, C, M;
    this->ComputeImplicitJacobian (t, u, v, a, K, C, M);

    SG::Algebra::Vec_t F;
    this->ComputeExplicitValue (t, u, v, F);

    // 计算右侧项 R = F - K * u - C * v
    SG::Algebra::Vec_t R = F - (K * u + C * v);

    // 求解加速度 a = M^(-1) * R
    SG::Algebra::TPardiso<SG::Algebra::Real_t> solver;
    solver.compute (M);
    solver.solve (R, _OUT a);
    // 1.lump 2. consistenct
}

// 在 MModel 中实现 ComputeImplicitValue 和 ComputeExplicitValue 方法
void MModel::ComputeImplicitValue (const SG::Algebra::Real_t& t,
                                   const SG::Algebra::Vec_t&  u,
                                   const SG::Algebra::Vec_t&  v,
                                   const SG::Algebra::Vec_t&  a,
                                   SG::Algebra::Vec_t&        f)
{}

void MModel::ComputeFValue (const SG::DataStructure::Common::Id subcaseId)
{
    for (auto&& dynaLoad : m_timeDynaLoadSet)
    {
        dynaLoad->computeA (m_dbServiceFactory, subcaseId);
    }

    for (auto&& dynaLoad : m_freqDynaLoadSet)
    {
        dynaLoad->computeA (m_dbServiceFactory, subcaseId);
    }
}

void MModel::ComputeExplicitValue (const SG::Algebra::Real_t& t,
                                   const SG::Algebra::Vec_t&  u,
                                   const SG::Algebra::Vec_t&  v,
                                   SG::Algebra::Vec_t& _OUT   g)
{
    auto&       assembly = SG::FEMSolver::AssemblyEigen::GetInstance ();
    std::size_t n        = assembly.GetFreeDOFSize ();
    g                    = SG::Algebra::Vec_t (n);
    g.setZero ();

    for (auto&& dynaLoad : m_timeDynaLoadSet)
    {
        g += dynaLoad->Interplate (m_dbServiceFactory, t);
    }
}

void MModel::ComputeImplicitJacobian (const SG::Algebra::Real_t& t,
                                      const SG::Algebra::Vec_t&  u,
                                      const SG::Algebra::Vec_t&  v,
                                      const SG::Algebra::Vec_t&  a,
                                      SG::Algebra::SpMat_t&      dfdu,
                                      SG::Algebra::SpMat_t&      dfdv,
                                      SG::Algebra::SpMat_t&      dfda)
{
    this->ComputeStiffnessMatrix (_OUT dfdu);  // K
    this->ComputeDampingMatrix (_OUT dfdv);    // C
    this->ComputeMassMatrix (_OUT dfda);       // M
}

// void MModel::ComputeExplicitJacobian (const SG::Algebra::Real_t& t,
//                                       const SG::Algebra::Vec_t&  u,
//                                       const SG::Algebra::Vec_t&  v,
//                                       SG::Algebra::SpMat_t&      dgdu,
//                                       SG::Algebra::SpMat_t&      dgdv)
// {}

void MModel::ComputeEigenValue (int ModesOrder, SG::Algebra::Vec_t& eigen, SG::Algebra::Mat_t& phi)
{
    SG::Algebra::SpMat_t K, M;
    this->ComputeStiffnessMatrix (K);
    this->ComputeMassMatrix (M);

    // 求解特征值和特征向量
    SG::Algebra::TArpack<SG::Algebra::Real_t> arpack;

    auto opt  = arpack.getDefaultOption ();
    opt.nev   = ModesOrder;
    opt.which = "LM";
    arpack.setOption (opt);

    std::vector<SG::Algebra::Real_t> eigenvalues (ModesOrder);
    std::vector<SG::Algebra::Vec_t>  eigenvectors (ModesOrder);

    // ToDO arpack 目前采用子空间迭代, 无法计算全部特征值
    if (K.cols () > ModesOrder)
    {
        arpack.solve (K, M, eigenvalues, eigenvectors);
    }
    else
    {
        const auto      n = K.rows ();
        Eigen::MatrixXd K1 (n, n);

        Eigen::MatrixXd M1 (n, n);
        for (int r = 0; r < n; r++)
        {
            for (int l = 0; l < n; l++)
            {
                const auto Kij = K.coeff (r, l);
                K1 (r, l)      = Kij;
                K1 (l, r)      = Kij;

                const auto Mij = M.coeff (r, l);
                M1 (r, l)      = Mij;
                M1 (l, r)      = Mij;
            }
        }

        // 使用eigen 的特征值计算器计算小规模问题(求解全部特征值)
        Eigen::GeneralizedEigenSolver<Eigen::MatrixXd> eigSolver;
        // std::cout << "K1 = \n" << K1 << std::endl;
        // std::cout << "M1 = \n" << M1 << std::endl;
        eigSolver.compute (K1, M1);
        const auto D_origin = eigSolver.eigenvalues ();
        const auto V_origin = eigSolver.eigenvectors ();
        // std::cout << "D_origin = \n" << D_origin << std::endl;
        // std::cout << "V_origin = \n" << V_origin << std::endl;
        // 创建索引数组用于排序
        std::vector<int> indices (D_origin.size ());
        std::iota (indices.begin (), indices.end (), 0);

        // 根据特征值的实部从小到大排序索引
        std::sort (indices.begin (), indices.end (), [&D_origin] (int i, int j) { return D_origin (i).real () < D_origin (j).real (); });

        // 根据排序后的索引重新排列特征值和特征向量
        Eigen::VectorXcd D (D_origin.size ());
        Eigen::MatrixXcd V (V_origin.rows (), V_origin.cols ());

        for (int i = 0; i < indices.size (); ++i)
        {
            D (i)     = D_origin (indices[i]);
            V.col (i) = -V_origin.col (indices[i]);
        }
        // 现在 D_sorted 和 V_sorted 就是按特征值从小到大排列的结果
        // std::cout << "Sorted eigenvalues:\n" << D << std::endl;
        // std::cout << "Sorted eigenvector:\n" << V << std::endl;
        for (int i = 0; i < V.cols (); ++i)
        {
            std::complex<double> norm_factor =
                std::sqrt ((V.col (i).adjoint () * M1.cast<std::complex<double>> () * V.col (i)) (0, 0));
            V.col (i) /= norm_factor;
        }
        // std::cout << "质量归一化后的特征向量 = \n" << V<< std::endl;
        for (int i = 0; i < opt.nev; ++i)
        {
            eigenvalues[i]  = D[i].real ();
            eigenvectors[i] = SG::Algebra::Vec_t (n);
            for (int j = 0; j < n; j++)
            {
                eigenvectors[i][j] = V (j, i).real ();
            }
        }
    }

    // 计算广义质量和广义刚度
    SG::DataStructure::Common::Real gMass, gStiff;
    eigen.resize(opt.nev);
    phi.resize(K.rows(), opt.nev);
    for (int i = 0; i < opt.nev; ++i)
    {
        eigen (i) = eigenvalues[i];

        const auto& eigV = eigenvectors[i];
        phi.col (i)      = eigV;
        const auto eigVT = eigV.transpose ();
        gMass            = eigVT * M * eigV;
        gStiff           = eigVT * K * eigV;
        this->SetEigenValue (i + 1, eigen (i), gMass, gStiff, 0, 0);
    }
}

Status MModel::QueryFrenqencePoint (Id tbId, std::vector<SG::DataStructure::Common::Real>& _OUT freSet)
{
    std::vector<SG::DataStructure::FEM::FREQ1Data> FREQ1Data;
    m_dbServiceFactory->get<IFREQ1Service> ()->FindBySId (tbId, FREQ1Data);
    freSet = SG::DataStructure::FEM::getfreqVector (FREQ1Data);
    if (freSet.empty ())
    {
        SGError << "can't find FREQ1Data  \n";
        return Status::STATUS_Failure;
    }

    return Status::STATUS_Success;
}

void MModel::ComputeDynaStiffness (SG::Algebra::ComplexSpMat_t& _OUT Kdyna)
{
    SG::Algebra::SpMat_t K;
    this->ComputeStiffnessMatrix (K);
    SG::Utility::SGParameter paraG = m_dbServiceFactory->get<IParameterService> ()->GetSolverParameter ();
    SG::Algebra::Real_t      G{ 0.0 };
    if (paraG.Has ("G"))
    {
        G = std::stof (paraG["G"].GetString ());
    }
    Kdyna = (Complex_t (1.0, G) * K);
}

void MModel::ComputeFreqForce (const SG::Algebra::Real_t&      t,
                               const SG::Algebra::Vec_t&       u,
                               const SG::Algebra::Vec_t&       v,
                               SG::Algebra::ComplexVec_t& _OUT Pf)
{
    auto&       assembly = SG::FEMSolver::AssemblyEigen::GetInstance ();
    std::size_t n        = assembly.GetFreeDOFSize ();
    Pf                   = SG::Algebra::ComplexVec_t (n);
    Pf.setZero ();

    for (auto&& load : m_freqDynaLoadSet)
    {
        Pf += load->Interplate (m_dbServiceFactory, t);
    }
}

void MModel::GetNodalSolution (const Variable& variable, SG::Algebra::Vec_t& value, const TimeIndex& t)
{
    SGWarn << __FILE__ << ":" << __LINE__ << " NOT IMPLEMENTED";
}

void MModel::SetNodalSolution (const Variable& variable, SG::Algebra::Vec_t& value, const TimeIndex& t)
{
    SGWarn << __FILE__ << ":" << __LINE__ << " NOT IMPLEMENTED";
}

void MModel::GetNodalSolution (const Variable& variable, SG::Algebra::Mat_t& data, const TimeIndex& t)
{
    SGWarn << __FILE__ << ":" << __LINE__ << " NOT IMPLEMENTED";
}

void MModel::SetNodalSolution (const Variable& variable, SG::Algebra::Mat_t& data, const TimeIndex& t)
{
    SGWarn << __FILE__ << ":" << __LINE__ << " NOT IMPLEMENTED";
}

void MModel::GetDisplacement (const SG::DataStructure::Common::Id subcaseId, const double t, SG::Algebra::Vec_t& V) const
{
    auto dataService = m_dbServiceFactory->get<IDISPLACEMENTService> ();

    std::vector<DISPLACEMENTData> nastranFormatData;
    if (Status::STATUS_Failure == dataService->Find ({ m_analysisType, subcaseId, t }, nastranFormatData))
    {
        return;
    }

    // 3. 数据转换
    std::vector<SG::DataStructure::Common::Real> realVec;
    SG::FEMSolver::DataAdapter::UnpackGridField (nastranFormatData, _OUT realVec);

    V.resize (realVec.size ());
    // 将realVec的数据赋值给V
    for (size_t i = 0; i < realVec.size (); ++i)
    {
        V[i] = realVec[i];
    }
}

// 目前仅仅考虑所有Free自由度的数值
SG::DataStructure::Common::Status
MModel::GetFreaRespMatrix (SG::DataStructure::Common::Real f, const IdArray& subcaseIds, SG::Algebra::ComplexMat_t& M)
{
    const auto k = subcaseIds.size ();
    if (0 == k)
    {
        return Status::STATUS_Failure;
    }

    auto&      assembly  = SG::FEMSolver::AssemblyEigen::GetInstance ();
    const auto nFreeDof = assembly.GetFreeDOFSize ();
    M                    = SG::Algebra::ComplexMat_t (nFreeDof, k);
    M.setZero ();

    SG::Algebra::ComplexVec_t Hk;
    std::vector<int> nodeGlobalDofIds (6);
    static const ::SG::DataStructure::FEM::DofSet DOFS{
    SG::DataStructure::FEM::DOFDirections::DOF_Ux,   SG::DataStructure::FEM::DOFDirections::DOF_Uy,
    SG::DataStructure::FEM::DOFDirections::DOF_Uz,   SG::DataStructure::FEM::DOFDirections::DOF_Rotx,
    SG::DataStructure::FEM::DOFDirections::DOF_Roty, SG::DataStructure::FEM::DOFDirections::DOF_Rotz
    };
    std::vector<DISPLACEMENT_CPLXData> nastranFormatData;
    auto dataService = m_dbServiceFactory->get<IDISPLACEMENT_CPLXService> ();
    for (size_t i = 0; i < k; ++i)
    {
        // GetDisplacement (subcaseIds[i], f, _OUT Hk);

        if (Status::STATUS_Failure == dataService->Find ({ m_analysisType, subcaseIds[i], f}, nastranFormatData))
        {
            return Status::STATUS_Failure;
        }

        // 3. 数据转换
        std::vector<std::complex<double>> realVec;
        SG::FEMSolver::DataAdapter::UnpackGridFieldCplx (nastranFormatData, _OUT realVec);
        // 仅仅考虑激活自由度
        Hk.resize (nFreeDof);
        // 将realVec的数据赋值给V
        for (const auto& nodeData : nastranFormatData)
        {
            assembly.ComputeGlobalDOFId (nodeData.m_id, DOFS, _OUT nodeGlobalDofIds);

            if (nodeGlobalDofIds[0] > -1 && nodeGlobalDofIds[0] < nFreeDof)
            {
                Hk[nodeGlobalDofIds[0]] = Complex_t(nodeData.m_XR, nodeData.m_XI);
            }
            if (nodeGlobalDofIds[1] > -1 && nodeGlobalDofIds[1] < nFreeDof)
            {
                Hk[nodeGlobalDofIds[1]] = Complex_t(nodeData.m_YR, nodeData.m_YI);
            }
            if (nodeGlobalDofIds[2] > -1 && nodeGlobalDofIds[2] < nFreeDof)
            {
                Hk[nodeGlobalDofIds[2]] = Complex_t(nodeData.m_ZR, nodeData.m_ZI);
            }
            if (nodeGlobalDofIds[3] > -1 && nodeGlobalDofIds[3] < nFreeDof)
            {
                Hk[nodeGlobalDofIds[3]] = Complex_t(nodeData.m_RXR, nodeData.m_RXI);
            }
            if (nodeGlobalDofIds[4] > -1 && nodeGlobalDofIds[4] < nFreeDof)
            {
                Hk[nodeGlobalDofIds[4]] = Complex_t(nodeData.m_RYR, nodeData.m_RYI);
            }
            if (nodeGlobalDofIds[5] > -1 && nodeGlobalDofIds[5] < nFreeDof)
            {
                Hk[nodeGlobalDofIds[5]] = Complex_t(nodeData.m_RZR, nodeData.m_RZI);
            }
        }


        if (static_cast<int> (Hk.size ()) != M.rows ())
        {
            SGError << "GetFreaRespMatrix: row mismatch";
            return Status::STATUS_Failure;
        }
        std::cout << "频响矩阵 = \n" << Hk << std::endl;
        M.col (static_cast<int> (i)) = Hk;
    }

    return Status::STATUS_Success;
}

void MModel::SetDisplacement (SG::Algebra::Vec_t V, const SG::DataStructure::Common::Id subcaseId, const double t)
{
    auto dataService = m_dbServiceFactory->get<IDISPLACEMENTService> ();

    // 位移恢复
    std::vector<SG::DataStructure::Common::Real> realVec;
    SG::FEMSolver::AssemblyEigen::GetInstance ().RecoveryDisplacement (V, realVec);

    // 转换数据格式
    std::vector<DISPLACEMENTData> nastranFormatData;
    SG::FEMSolver::DataAdapter::PackGridField (realVec, _OUT nastranFormatData);

    // 添加数据
    if (Status::STATUS_Failure == dataService->Add ({ m_analysisType, subcaseId, t }, nastranFormatData))
    {
        SGWarn << __FILE__ << ":" << __LINE__ << "Save " << typeid (EIGENVECTORData).name () << " is failed！";
    }
}

void MModel::SetDisplacement (const SG::Algebra::Vec_t& V, Id subcaseId, SG::Algebra::Real_t t, RandomResultType random)
{
    auto dataService = m_dbServiceFactory->get<IDISPLACEMENTService> ();

    // 位移恢复
    std::vector<SG::Algebra::Real_t> realVec;
    SG::FEMSolver::AssemblyEigen::GetInstance ().RecoveryDisplacement (const_cast<SG::Algebra::Vec_t&>(V), realVec);

    // // 转换数据格式
    std::vector<DISPLACEMENTData> nastranFormatData;
    SG::FEMSolver::DataAdapter::PackGridField (realVec, _OUT nastranFormatData);

    // // 添加数据
    if (Status::STATUS_Failure == dataService->Add ({ subcaseId, t, random }, nastranFormatData))
    {
        SGWarn << __FILE__ << ":" << __LINE__ << "Save " << typeid (EIGENVECTORData).name () << " is failed！";
    }
}

void MModel::GetDisplacement (const SG::DataStructure::Common::Id subcaseId, const SG::Algebra::Real_t t, SG::Algebra::ComplexVec_t& V) const
{
    auto dataService = m_dbServiceFactory->get<IDISPLACEMENT_CPLXService> ();

    std::vector<DISPLACEMENT_CPLXData> nastranFormatData;
    if (Status::STATUS_Failure == dataService->Find ({ m_analysisType, subcaseId, t }, nastranFormatData))
    {
        return;
    }

    // 3. 数据转换
    std::vector<std::complex<double>> realVec;
    SG::FEMSolver::DataAdapter::UnpackGridFieldCplx (nastranFormatData, _OUT realVec);

    V.resize (realVec.size ());
    // 将realVec的数据赋值给V
    for (size_t i = 0; i < realVec.size (); ++i)
    {
        V[i] = realVec[i];
    }
}

void MModel::GetVelocity (const SG::DataStructure::Common::Id subcaseId, const double t, SG::Algebra::Vec_t& V) const
{
    auto dataService = m_dbServiceFactory->get<IVELOCITYService> ();

    std::vector<DISPLACEMENTData> nastranFormatData;
    if (Status::STATUS_Failure == dataService->Find ({ m_analysisType, subcaseId, t }, nastranFormatData))
    {
        return;
    }

    // 3. 数据转换
    std::vector<SG::DataStructure::Common::Real> realVec;
    SG::FEMSolver::DataAdapter::UnpackGridField (nastranFormatData, _OUT realVec);

    V.resize (realVec.size ());
    // 将realVec的数据赋值给V
    for (size_t i = 0; i < realVec.size (); ++i)
    {
        V[i] = realVec[i];
    }
}

void MModel::GetAcceleration (const SG::DataStructure::Common::Id subcaseId, const double t, SG::Algebra::Vec_t& V) const
{
    auto dataService = m_dbServiceFactory->get<IACCELERATIONService> ();

    std::vector<DISPLACEMENTData> nastranFormatData;
    if (Status::STATUS_Failure == dataService->Find ({ m_analysisType, subcaseId, t }, nastranFormatData))
    {
        return;
    }

    // 3. 数据转换
    std::vector<SG::DataStructure::Common::Real> realVec;
    SG::FEMSolver::DataAdapter::UnpackGridField (nastranFormatData, _OUT realVec);

    V.resize (realVec.size ());
    // 将realVec的数据赋值给V
    for (size_t i = 0; i < realVec.size (); ++i)
    {
        V[i] = realVec[i];
    }
}

void MModel::SetDisplacement (SG::Algebra::ComplexVec_t& V, const SG::DataStructure::Common::Id subcaseId, const double t)
{
    auto dataService = m_dbServiceFactory->get<IDISPLACEMENT_CPLXService> ();

    // 位移恢复
    std::vector<std::complex<double>> realVec;
    SG::FEMSolver::AssemblyEigen::GetInstance ().RecoveryDisplacement (V, realVec);

    // 转换数据格式
    std::vector<DISPLACEMENT_CPLXData> nastranFormatData;
    SG::FEMSolver::DataAdapter::PackGridFieldCplx (realVec, _OUT nastranFormatData);

    // 添加数据
    if (Status::STATUS_Failure == dataService->Add ({ this->m_analysisType, subcaseId, t }, nastranFormatData))
    {
        SGWarn << __FILE__ << ":" << __LINE__ << "Save " << typeid (EIGENVECTORData).name () << " is failed！";
    }
}

void MModel::SetVelocity (SG::Algebra::Vec_t V, const SG::DataStructure::Common::Id subcaseId, const double t)
{
    auto dataService = m_dbServiceFactory->get<IVELOCITYService> ();

    // 位移恢复
    std::vector<SG::DataStructure::Common::Real> realVec;
    SG::FEMSolver::AssemblyEigen::GetInstance ().RecoveryDisplacement (V, realVec);

    // 转换数据格式
    std::vector<DISPLACEMENTData> nastranFormatData;
    SG::FEMSolver::DataAdapter::PackGridField (realVec, _OUT nastranFormatData);

    // 添加数据
    if (Status::STATUS_Failure == dataService->Add ({ this->m_analysisType, subcaseId, t }, nastranFormatData))
    {
        SGWarn << __FILE__ << ":" << __LINE__ << "Save " << typeid (EIGENVECTORData).name () << " is failed！";
    }
}

void MModel::SetAcceleration (SG::Algebra::Vec_t V, const SG::DataStructure::Common::Id subcaseId, const double t)
{
    auto dataService = m_dbServiceFactory->get<IACCELERATIONService> ();

    // 位移恢复
    std::vector<SG::DataStructure::Common::Real> realVec;
    SG::FEMSolver::AssemblyEigen::GetInstance ().RecoveryDisplacement (V, realVec);

    // 转换数据格式
    std::vector<DISPLACEMENTData> nastranFormatData;
    SG::FEMSolver::DataAdapter::PackGridField (realVec, _OUT nastranFormatData);

    // 添加数据
    if (Status::STATUS_Failure == dataService->Add ({ this->m_analysisType, subcaseId, t }, nastranFormatData))
    {
        SGWarn << __FILE__ << ":" << __LINE__ << "Save " << typeid (EIGENVECTORData).name () << " is failed！";
    }
}

void MModel::GetVelocity (SG::Algebra::Mat_t& V, const SG::DataStructure::Common::Id subcaseId, const TimeIndex& t) const
{
    SGWarn << __FILE__ << ":" << __LINE__ << " NOT IMPLEMENTED";
}

void MModel::SetVelocity (SG::Algebra::Mat_t& V, const SG::DataStructure::Common::Id subcaseId, const TimeIndex& t)
{
    SGWarn << __FILE__ << ":" << __LINE__ << " NOT IMPLEMENTED";
}

void MModel::GetAcceleration (SG::Algebra::Mat_t& A, const SG::DataStructure::Common::Id subcaseId, const TimeIndex& t) const
{
    SGWarn << __FILE__ << ":" << __LINE__ << " NOT IMPLEMENTED";
}

void MModel::SetAcceleration (SG::Algebra::Mat_t& A, const SG::DataStructure::Common::Id subcaseId, const TimeIndex& t)
{
    SGWarn << __FILE__ << ":" << __LINE__ << " NOT IMPLEMENTED";
}

void Kernel::MModel::GetEigenVector (std::vector<SG::DataStructure::Common::Real>& V,
                                     const SG::DataStructure::Common::Id           subcaseId,
                                     const TimeIndex&                              t) const
{
    auto dataService = m_dbServiceFactory->get<IEIGENVECTORService> ();

    std::vector<EIGENVECTORData> nastranFormatData;
    if (Status::STATUS_Failure == dataService->Find (m_analysisType, subcaseId, t, nastranFormatData))
    {
        return;
    }

    // 3. 数据转换
    SG::FEMSolver::DataAdapter::UnpackGridField (nastranFormatData, _OUT V);
}

void Kernel::MModel::SetCrossPSD (const std::vector<SG::DataStructure::Post::CROSS_PSDData>& datas,
                                  const Id                                                   randomSid,
                                  const SG::Algebra::Real_t                                  freq)
{
    auto dataService = m_dbServiceFactory->get<ICROSS_PSDService> ();

    // 添加数据
    if (Status::STATUS_Failure == dataService->Add ({ randomSid, freq, RandomResultType::PSDS }, const_cast<std::vector<CROSS_PSDData>&> (datas)))
    {
        SGWarn << __FILE__ << ":" << __LINE__ << "Save " << typeid (CROSS_PSDData).name () << " is failed！";
    }
}


void Kernel::MModel::SetEigenVector (SG::Algebra::Mat_t&                 eigenVectorM,
                                     SG::Algebra::Vec_t                  eigenValue,
                                     const SG::DataStructure::Common::Id subcaseId)
{
    std::vector<SG::DataStructure::Common::Real> eigenVector;

    auto dataService = m_dbServiceFactory->get<IEIGENVECTORService> ();

    // 结果
    std::vector<EIGENVECTORData> nastranFormatData;

    for (auto c = 0; c < eigenVectorM.cols (); c++)
    {
        // 阶
        int  mode = c + 1;
        Real freo = eigenValue (c);

        // 位移恢复
        SG::FEMSolver::EigenVector phi = eigenVectorM.col (c);
        SG::FEMSolver::AssemblyEigen::GetInstance ().RecoveryDisplacement (phi, eigenVector);

        // 转换数据格式
        SG::FEMSolver::DataAdapter::PackGridField (eigenVector, _OUT nastranFormatData);

        // 添加数据
        if (Status::STATUS_Failure == dataService->Add (this->m_analysisType, subcaseId, mode, freo, nastranFormatData))
        {
            SGWarn << __FILE__ << ":" << __LINE__ << "Save " << typeid (EIGENVECTORData).name () << " is failed！";
        }
        nastranFormatData.clear ();
    }
}

void Kernel::MModel::GetEigenVector (const Id subcaseId, const TimeIndex& t, _OUT std::vector<SG::Algebra::Complex_t>& V) const
{
    auto dataService = m_dbServiceFactory->get<IEIGENVECTOR_CPLXService> ();

    std::vector<EIGENVECTOR_CPLXData> nastranFormatData;
    if (Status::STATUS_Failure == dataService->Find (m_analysisType, subcaseId, t, nastranFormatData))
    {
        return;
    }

    // 3. 数据转换
    SG::FEMSolver::DataAdapter::UnpackGridFieldCplx (nastranFormatData, _OUT V);
}

void Kernel::MModel::SetEigenVector (SG::Algebra::ComplexMat_t& complexMat, SG::Algebra::ComplexVec_t eigenValue, const Id subcaseId)
{
    std::vector<std::complex<double>> eigenVector;
    auto                              dataService = m_dbServiceFactory->get<IEIGENVECTOR_CPLXService> ();

    // 结果
    std::vector<EIGENVECTOR_CPLXData> nastranFormatData;

    for (auto c = 0; c < complexMat.cols (); c++)
    {
        // 阶
        int  mode = c + 1;
        auto freo = eigenValue (c);

        // 位移恢复
        SG::Algebra::ComplexVec_t phi = complexMat.col (c);
        SG::FEMSolver::AssemblyEigen::GetInstance ().RecoveryDisplacement (phi, eigenVector);

        // 转换数据格式
        SG::FEMSolver::DataAdapter::PackGridFieldCplx (eigenVector, _OUT nastranFormatData);

        // 添加数据
        if (Status::STATUS_Failure == dataService->Add (m_analysisType, subcaseId, mode, freo, nastranFormatData))
        {
            SGWarn << __FILE__ << ":" << __LINE__ << "Save " << typeid (EIGENVECTOR_CPLXData).name () << " is failed！";
        }
        nastranFormatData.clear ();
    }
}

void Kernel::MModel::SetEigenValue (const int t, const Real eigen, const Real mass, const Real stiff, const I64 resFlag, const I64 fldFlag)
{
    auto dataService = m_dbServiceFactory->get<IEIGENVALUEService> ();

    SG::DataStructure::Post::EIGENVALUEData data;
    data.m_ORDER     = t;
    data.m_MODE      = t;
    data.m_EIGEN     = eigen;
    const auto omega = std::sqrt (fabs (eigen));
    data.m_OMEGA     = omega;
    data.m_FREQ      = omega / (2.0 * M_PI);
    data.m_MASS      = mass;
    data.m_STIFF     = stiff;
    data.m_RESFLG    = resFlag;
    data.m_FLDFLG    = fldFlag;

    if (Status::STATUS_Failure == dataService->Add (data))
    {
        SGWarn << __FILE__ << ":" << __LINE__ << " Save EIGENVALUEData is failed！";
    }
}

void Kernel::MModel::GetEigenValue (const I64  subcaseId,
                                    const int  t,
                                    _OUT Real& eigen,
                                    _OUT Real& omega,
                                    _OUT Real& freq,
                                    _OUT Real& mass,
                                    _OUT Real& stiff,
                                    _OUT I64&  resFlag,
                                    _OUT I64&  fldFlag) const
{
    auto dataService = m_dbServiceFactory->get<IEIGENVALUEService> ();

    SG::DataStructure::Post::EIGENVALUEData data;
    if (Status::STATUS_Success == dataService->Find (subcaseId, t, _OUT data))
    {
        // 数据赋值
        eigen   = data.m_EIGEN;
        omega   = data.m_OMEGA;
        freq    = data.m_FREQ;
        mass    = data.m_MASS;
        stiff   = data.m_STIFF;
        resFlag = data.m_RESFLG;
        fldFlag = data.m_FLDFLG;
    }
    else
    {
        SGWarn << __FILE__ << ":" << __LINE__ << " Get EIGENVALUEData is failed！";
    }
}

void Kernel::MModel::SetEigenValue (const int t, const SG::Algebra::Complex_t eigen, const Real freq, const Real damp, const Real spinSpeed)
{
    auto dataService = m_dbServiceFactory->get<IEIGENVALUE_CPLXService> ();

    SG::DataStructure::Post::EIGENVALUE_CPLXData data;
    data.m_ORDER      = t;
    data.m_MODE       = t;
    data.m_REIGEN     = eigen.real ();
    data.m_IEIGEN     = eigen.imag ();
    data.m_FREQ       = freq;
    data.m_DAMP       = damp;
    data.m_SPIN_SPEED = spinSpeed;

    if (Status::STATUS_Failure == dataService->Add (data))
    {
        SGWarn << __FILE__ << ":" << __LINE__ << "Save EIGENVALUE_CPLXData is failed！";
    }
}

void Kernel::MModel::GetEigenValue (
    const I64 subcaseId, const int t, _OUT SG::Algebra::Complex_t& eigen, _OUT Real& freq, _OUT Real& damp, _OUT Real& spinSpeed) const
{
    auto dataService = m_dbServiceFactory->get<IEIGENVALUE_CPLXService> ();

    SG::DataStructure::Post::EIGENVALUE_CPLXData data;
    if (Status::STATUS_Success == dataService->Find (subcaseId, t, _OUT data))
    {
        // 数据赋值
        eigen     = { data.m_REIGEN, data.m_IEIGEN };
        freq      = data.m_FREQ;
        damp      = data.m_DAMP;
        spinSpeed = data.m_SPIN_SPEED;
    }
}

void MModel::ComputeOutputData (const TimeIndex& t)
{
    SGWarn << __FILE__ << ":" << __LINE__ << " NOT IMPLEMENTED";
}

void MModel::Finalize ()
{
    SGWarn << __FILE__ << ":" << __LINE__ << " NOT IMPLEMENTED";
}
